Helicoidal instability of a scroll vortex in three-dimensional reaction-diffusion systems 
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We study the dynamics of scroll vortices in excitable reaction-diffusion systems analytically and 
numerically. We demonstrate that intrinsic three-dimensional instability of a straight scroll leads to 
the formation of helicoidal structures. This behavior originates from the competition between the 
scroll curvature and unstable core dynamics. We show that the obtained instability persists even 
beyond the meander core instability of two-dimensional spiral wave. 
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Spiral waves arising in diverse physical, chemical, and 
biological systems are now one of the paradigms of non- 
equilibrium dynamical phenomena [[l] . Examples include 
the Belousov- Zhabotinsky (BZ) reaction jl and catalytic 
oxidation of CO on Pt substrates || , concentration waves 
in colonies of aggregating amoebae 0], waves in cardiac 
tissue |^], and many others. These seemingly unrelated 
phenomena share a common feature: they often allow 
for a description within the framework of two-component 
reaction-diffusion type systems, for which spiral solutions 
are generic. 

The three-dimensional analog of a spiral wave is a scroll 
vortex, which can be represented by translating the spi- 
ral wave along the third direction. Thus, the point sin- 
gularity of two-dimensional spiral wave (tip) develops 
into a line singularity (vortex filament). The dynam- 
ics of the vortex filaments in reaction diffusion systems 
has attracted a great deal of attention in connection 
with sudden heart fibrillation, where it is believed that 
scroll and ring vortices play a crucial role [^]|| . Extensive 
numerical simulations [^|-^| and recent experiments on 
gel-immobilized BZ reaction show that in a wide 

range of parameters the scrolls are unstable and may as- 
sume helicoidal or even more complicated dynamic con- 
figurations. However, the theoretical analysis performed 
first by J. Keener, predicts an ultimate collapse of vortex 
rings and stability of straight vortex filaments jl2| . This 
conclusion was drawn on the basis of multiscale analysis, 
which shows that a vortex ring is similar to an elastic line 
with a positive line tension. The persistence of nontrivial 
vortex configurations and turbulence in reaction-diffusion 
systems was attributed to a negative line tension of the 
filament ^]. 

In this Letter we demonstrate, on the basis of numer- 
ical and analytical calculations, that the formation of 
helicoidal vortices can be related to the intrinsic three- 
dimensional instability of a straight scroll, caused by a 
nontrivial response of the filament core to a bending of 
the filament. We show that the limit of this instabil- 
ity, resulting in the formation of spontaneous helicoidal 



vortices, goes beyond the corresponding two-dimensional 
core meander instability. 

The dynamics of a scroll vortex can be consistently de- 
scribed by the two-component reaction-diffusion system 



dtu = eV 2 u 



f(u,v) 



d t v — 5eV 2 v + g(u, v) , 



(1) 
(2) 



where e is a small positive parameter, and 5 = D v /D u 
is the ratio of diffusion coefficients of the variables v and 
u. The functions f(u, v) and g(u, v) are chosen so as 
to make Eqs. ([T|)-(p|) excitable [Q. In a wide range of 
parameters Eqs. ([j])-(||) have a spiral wave solution in 
two dimensions and a scroll vortex in three dimensions. 

Let us derive the equation of motion for the core of 
a weakly-curved scroll vortex subject to the meander- 
ing instability [^|,^3). In two dimensions this instability 
was studied both via numerical simulations of the system 
(Q)-(||) , and using the numerical solution of the linearized 
problem |l4|,[l5|]. Recently it was proven that the spiral 
interface undergoes a core-meander instability via a su- 
percritical Hopf bifurcation, as the diffusion coefficient 
of the slow field S decreases |l6|] . The mechanism of the 
meandering instability in a certain parameter limit has 
been elaborated on in a recent work |l7|l . 

The core meander was described by five equations for 
the frequency, coordinates, and the velocity of the spiral 
tip | fl5l| . These equations can be considerably simplified 
if, instead of the spiral tip, one considers the instant cen- 
ter of spiral rotation. At the threshold of instability one 
obtains a single complex Landau-type equation for the 
"complex" velocity of the rotation center C = c x + ic y , 
where c x ^ v are the components of the velocity: 



d t C 



yC-/3\C\ 2 C . 



(3) 



Here a and (3 are complex coefficients, a — ct\ + ia2 , /3 = 
Pi + ?'/?2 , that have to be determined numerically. For 
ct\ < the symmetry center is stable, and the spiral 
makes a pure rotation. When ol\ becomes positive, the 
fixed point solution of Eq. (0) loses stability, and the 
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rotation center itself performs a circular motion, which 
implies the meandering (composite rotation) of the spiral 
tip ]l8| . As follows from Eq. (||), the absolute value of 
the velocity of the rotation center, in the saturated mean- 
dering regime, is Cq = y/ai/\/3i\, and the corresponding 
rotation frequency is ujq = a-i — U\f$ilfi\ ■ 

Consider now a three-dimensional weakly curved scroll 
(the curvature k in the third dimension is much smaller 
then the local curvature of the spiral front). In this limit 
one expands the Laplacian in equation (|l|) as follows, 
V 2 w V?, D — kN ■ V, where V| D is the Laplacian in 
the cross-section, and N is the unit vector pointing to- 
ward the center of the filament curvature. Therefore, the 
curvature in the third dimension plays the role of an ad- 
vective field directed along N and causing the drift of 
the filament [fil|. To lowest order, the curvature k enters 
into the equation of motion (^) linearly: 



d t C = aC 



(4) 



where 7 = 71 + 172 is a (complex) constant which can 
be determined numerically from two-dimensional simu- 
lations or analytically in a large core limit jL7|. 

Equation (^) readily implies the short-wavelength in- 
stability of a straight filament. Indeed, for an almost 
straight filament parallel to the z-axis, the curvature 
vector kN is simply (x zz ,y zz ). Using C — dtx , where 
x = x + iy, we obtain from Eq. (Q), for a periodic per- 
turbation of the filament x(z) ~ exp[A(fe)t + ikz] , that 
A 2 = a (A + jk 2 ) . The latter equation has the following 
unstable solution 



X(k) 



a 
2 



ajk 2 



(5) 



Eq. (jq) represents the eigenvalue of the linearized prob- 
lem (^) and is valid only for slightly curved filaments, 
i.e. for small wavenumbers k. If this near-threshold in- 
stability saturates, the saturated structure corresponds 
to the most unstable mode, which implies the formation 
of stable helicoidal vortices. 

By analogy with the vortex filaments in the CGLE [|o| , 
we can expect that, for k > , the growth rate of the ob- 
tained three-dimensional instability, Re[X(k)], substan- 
tially exceeds that of the two-dimensional meandering 
instability, i?e[A(0)] = Re[a). In that case the filament 
curvature plays a destabilizing role in the dynamics of the 
filament. To analyze the stability, one can find the coeffi- 
cients a and 7 from the dynamics of the two-dimensional 
system (|l])-(||), and then substitute them into Eq. (|J). 

To determine a and 7 numerically, we have simulated 
Eqs. (|l|)-(||) in two dimensions using the EZ-spiral code 
of D. Barkley Ej]. We have chosen Barkley's model 
with the functions f(u,v) = u(u — l)[u — u t h(v)] and 
g(u, v) = u — v, where u t h{v) — (v + b)/a. To deter- 
mine numerically a, (3 and 7, one should start with an 
unstable rigidly rotating spiral, which is not available in 



numerical simulations. To overcome this difficulty, we 
used the following approach. We started with a spiral 
with already developed meander. Then we applied a lo- 



calized control technique, developed in Ref. |22|, to turn 
the spiral motion to a pure rotation around its symmetry 
center. Following Ref. plj , we applied a pinning source 
to equation (^|) , in the form of a localized inhibiting term 
—fih(r — r*o) , where rg is the coordinate of the domain 
center. Here /1 plays the role of a control parameter and 
is governed by the equation 



d t [i = -aifi + bi[v(r a ) - v ] , 



(6) 



where vq is the expected value of the slow field at the spi- 
ral center, and a\, b\ are coefficients that should be prop- 
erly chosen. We have elaborated the controlling tech- 
nique by allowing Vq to vary adiabatically (slowly com- 
pared to jj) as follows 



d t v = a 2 ^ + b 2 [v(r ) - v ] 



(7) 



with another pair of coefficients 02, &2- In this way, we 
have obtained a purely rotating spiral with the parameter 
values corresponding to a meandering regime. The final 
value of the control parameter was as small as fj, = I0~ 4 , 
which made the source term added to equation (Q) a small 
perturbation. After "switching off" the control (setting 
fi = 0), we allow the linear meander instability to de- 
velop. 




Time 

FIG. 1. Coordinates of spiral rotation center as func- 
tions of time, in the regime of linear development of me- 
andering instability. Parameters of Barkley's model are 
a — 0.66 , b = 0.01 , e = 1/50 . System size is 10 x 10, number 
of grid-points is 81 x 81. Time step dt w 0.0016. 

In order to find the constant a, we fitted the trajectory 
of the spiral tip to the developing meander, in the form 
x(t) = xq + x\ sin(cjii + <p±) + xi exp(ait) sin(a2t + (fn) , 
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where xq is the average position of the spiral tip, oji is 
the main spiral rotation frequency, and a = oi\ + ia 2 . 
The same fitting was carried out for y(t). Thus, we have 
retrieved the motion of the center of spiral rotational 
symmetry. The coordinates of the center, x c and y c , as 
functions of time are given in Fig. [I]. We have found that, 
for the set of parameters given in the caption to Fig. [j], 
ati « 0.012 and a 2 ~ 0.336 . To find the constant 7 
we have applied a homogeneous advective field to equa- 
tion ([!]), along the x-axis, as we have done in Ref. Jl9| . 
Then 7 is determined from fitting the spiral tip trajectory 
to the meandering in the above form with an additional 
term jiEt (72-Et), in the expression for x(t) (y(t) ), due 
to the drift of the spiral. The result for 7 is 71 rs 1.312 
and 72 w 0.154. In principle, it is easy to obtain the 
value of the parameter (3 by matching the spiral center 
trajectory to the saturated meander. However, we do not 
do it because it is not important for the studied mecha- 
nism of the instability. 




time (the transient from random initial conditions is very 
long). 

As we can see from Fig. pi the k — » asymptotics of 
.Re [A] and 7m [A] obtained numerically are well matched 
by our theoretical prediction. As k increases, however, 
the small curvature approximation used for the deriva- 
tion of (^) no longer works, and the prediction becomes 
invalid. This is why we see the divergence of the theo- 
retical curves from the numerical ones for larger k. 




FIG. 2. Three-dimensional filament with the saturated in- 
stability. The system size is 10 x 10x40, number of grid-points 
is 81 x 81 x 320. Other parameters are the same as in Fig [!} 
The initial perturbation is of the 3rd harmonics, with the 
wavenumber k « 0.47 . 

We have performed numerical simulations of the three- 
dimensional Eqs. (Q)-(0). We have studied the behavior 
of an almost straight scroll. We have prepared a two- 
dimensional purely rotating (unstable) spiral, with the 
parameter values in the meandering regime, using the 
elaborated controlling technique described above, and 
translated it along the third dimension to build an ini- 
tial scroll. Then, a periodic perturbation ~ exp[ifc,z] has 
been applied. The three-dimensional plot of the helix, ob- 
tained as the result of the instability saturation, is given 
in Fig. ^. In Fig. || we have plotted the real and imag- 
inary parts of the eigenvalue A versus k, both for the 
theoretical prediction given by (0), and for the results 
of the simulations. As we have predicted, the observed 
instability of the filament with initial perturbation of a 
finite k appears to be substantially stronger than the two- 
dimensional core meander instability ( k — ) . We have 
also tried initial conditions other than purely periodic, 
and have always found a saturated helix. Almost peri- 
odic initial conditions have been mainly used to cut CPU 
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FIG. 3. The real (a) and imaginary (b) parts of the eigen- 
value X(k) of the three-dimensional instability obtained ana- 
lytically (solid line) and numerically (dashed line). The sys- 
tem parameters are the same as in Fig ^. 

We have also studied numerically the dynamics of fil- 
aments in the space of parameters a and b of Barkley's 
model We have performed the stability analysis of 
the harmonic with the growth rate near the maximal one 
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(see Fig. ||(a) ). In Fig. || we plot the bifurcation line of 
the obtained three-dimensional instability together with 
that of the core meander instability. We see that the cur- 
vature of the filament enhances its instability, so that it 
is unstable even in a domain of the a—b plane, for which 
the two-dimensional spiral is stable. 

In conclusion, we have demonstrated that the intrinsic 
three-dimensional instability of a straight scroll leads to 
the formation of stable helicoidal structures. We have 
shown that persistent curved vortex configurations are 
not necessarily related to a " negative line tension" of the 
filament, but originate from the underdamped core dy- 
namics. Our current analysis is restricted to untwisted 
scrolls. We expect that twisted scrolls are unstable even 
in a wider space of the parameters, and exhibit in general 
even more violent dynamics. Helicoidal structures with 
twist were observed recently in gel-immobilized BZ reac- 
tion jn]]. The twist was generated by an external tem- 
perature gradient along the scroll axis. The helicoidal 
instability was observed above specific value of the twist 
and disappeared when the temperature gradient (and, 
therefore, the twist) was removed. Presumably, by tun- 
ing the parameters of the BZ reaction, one may approach 
the limit of the intrinsic three-dimensional instability for 
untwisted scrolls. We can also speculate that the heli- 
coidal instability of scroll vortices in large aspect ratio 
reaction-diffusion systems is one of the mechanisms that 
supports more complex vortex configurations and drives 
turbulent states. 
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FIG. 4. Bifurcation lines of 2D (meander) instability 
(solid line) and 3D instability (dashed line) in the plane a-b of 
the Barkley's model parameters. The wavenumber k « 0.63 
of the initial perturbation corresponds to the 4th harmonics, 
near the maximal growth rate of the instability, according to 
Fig. H(a) . Other parameters are the same as in Fig. 0. 
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